Manuscript_CAPTAIN_analyses

Load libraries

library(here)
library(dplyr)
library(ggplot2)
library(sf)
library(rnaturalearth)
library(smoothr)
library(raster)
library(readr)
library(tidyr)
library(ggpubr)
library(betareg)
library(broom)
library(margins)
library(patchwork)
library(dplyr)
library(sf)
library(egg)
library(spatstat)
library(sf)
library(tidyverse)      # For data manipulation and visualization
library(rnaturalearthdata) # Additional natural earth data
library(biscale)        # For bivariate mapping
library(gridExtra)      # For arranging multiple plots
library(grid)           # For text elements in plots
library(colorspace)     # For color manipulation

Load the data

# Read the CAPTAIN2 EDGE2 RDS file
CAPTAIN2_EDGE2_data <- readRDS(here::here("Data/CAPTAIN2_EDGE_full_results_averaged_budget0.1_replicates50.rds"))

# Read the CAPTAIN2 FUSE RDS file
CAPTAIN2_FUSE_data <- readRDS(here::here("Data/CAPTAIN2_FUSE_res_full_results_averaged_budget0.1_replicates50.rds"))

# Read the CAPTAIN2 IUCN RDS file
CAPTAIN2_IUCN_data <- readRDS(here::here("Data/CAPTAIN2_IUCN_full_results_averaged_budget0.1_replicates50.rds"))

# Load one of your input raster files to extract the correct grid structure
# Use the same raster file for consistency
raster_file <- here::here("Data", "tif files continental", "Psammobatis_parvacauda.tif")

# Check if the file exists
if (!file.exists(raster_file)) {
  stop("Raster file not found. Please provide a valid path to one of your input raster files.")
}

# Load the raster
r <- raster(raster_file)

# Read the protected range fractions RDS file
protected_fractions <- readRDS(here::here("Data", "CAPTAIN2_protected_range_fractions_2ndrun.rds"))

# Read the continental shark conservation metrics CSV file
shark_metrics <- read_csv(here::here("Data", "continental_shark_conservation_metrics_10_harmonised_IUCN_categories.csv"))

# Load fishing data
load(here::here("Data", "Raw", "Predicted_Fishing_Hours_05Deg.Rdata"))

EDGE2 continental 0.1 budget

# Analyze non-zero cells in the prioritization output
cat("Analyzing cells with non-zero priority values in CAPTAIN2 EDGE2 output:\n")

Analyzing cells with non-zero priority values in CAPTAIN2 EDGE2 output:

total_cells <- nrow(CAPTAIN2_EDGE2_data)
nonzero_cells <- sum(CAPTAIN2_EDGE2_data$Priority > 0, na.rm = TRUE)
zero_cells <- sum(CAPTAIN2_EDGE2_data$Priority == 0, na.rm = TRUE)
na_cells <- sum(is.na(CAPTAIN2_EDGE2_data$Priority))

cat("Total cells in grid:", total_cells, "\n")

Total cells in grid: 232560 

cat("Cells with non-zero priority:", nonzero_cells, " (", round(nonzero_cells/total_cells*100, 2), "%)\n", sep="")

Cells with non-zero priority:5618 (2.42%)

cat("Cells with zero priority:", zero_cells, " (", round(zero_cells/total_cells*100, 2), "%)\n", sep="")

Cells with zero priority:226942 (97.58%)

cat("Cells with NA priority:", na_cells, " (", round(na_cells/total_cells*100, 2), "%)\n", sep="")

Cells with NA priority:0 (0%)

# Summary statistics of priority values
priority_summary <- summary(CAPTAIN2_EDGE2_data$Priority)
cat("\nSummary statistics of priority values:\n")

Summary statistics of priority values:

print(priority_summary)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
0.00000 0.00000 0.00000 0.02331 0.00000 1.00000 

# Distribution of non-zero priority values
nonzero_priority <- CAPTAIN2_EDGE2_data$Priority[CAPTAIN2_EDGE2_data$Priority > 0]
cat("\nDistribution of non-zero priority values:\n")

Distribution of non-zero priority values:

priority_quantiles <- quantile(nonzero_priority, probs = seq(0, 1, 0.1), na.rm = TRUE)
print(priority_quantiles)

  0%  10%  20%  30%  40%  50%  60%  70%  80%  90% 100% 
0.02 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 

# Load one of your input raster files to extract the correct grid structure
raster_file <- here::here("Data", "tif files continental", "Psammobatis_parvacauda.tif")

# Check if the file exists
if (!file.exists(raster_file)) {
  stop("Raster file not found. Please provide a valid path to one of your input raster files.")
}

# Load the raster
r <- raster(raster_file)

# Get the dimensions of the raster
nrows <- nrow(r)
ncols <- ncol(r)

# Confirm dimensions match expected values
if (nrows != 323 || ncols != 720) {
  warning("Raster dimensions don't match expected values. Proceeding with actual dimensions.")
}

# Create a grid of coordinates for each cell
coords <- as.data.frame(coordinates(r))
names(coords) <- c("Longitude", "Latitude")

# Add cell IDs (PUID) to the coordinates
coords$PUID <- 1:nrow(coords)

# Now join with the CAPTAIN2 data based on PUID
CAPTAIN2_EDGE2_data_with_coords <- CAPTAIN2_EDGE2_data %>%
  left_join(coords, by = "PUID")

# Check if the join worked correctly
if (sum(is.na(CAPTAIN2_EDGE2_data_with_coords$Longitude)) > 0) {
  warning("Some PUIDs from CAPTAIN2 data couldn't be matched to coordinates.")
}

# Filter to keep only cells with non-zero priority for faster plotting
CAPTAIN2_EDGE2_data_nonzero <- CAPTAIN2_EDGE2_data_with_coords %>%
  filter(Priority > 0) %>%
  filter(!is.na(Longitude), !is.na(Latitude))  # Remove any rows with missing coords

# Get world map data
world <- ne_countries(scale = "medium", returnclass = "sf")

# Define the McBryde-Thomas 2 projection
mcbryde_thomas_2 <- "+proj=mbt_s"

# Transform the dataset to sf object and project
CAPTAIN2_EDGE2_sf <- st_as_sf(
  CAPTAIN2_EDGE2_data_nonzero, 
  coords = c("Longitude", "Latitude"), 
  crs = crs(r, asText = TRUE)  # Use the raster's CRS
) %>%
  st_transform(crs = mcbryde_thomas_2)

# Project the world map
world_projected_CAPTAIN2_EDGE2 <- st_transform(world, crs = mcbryde_thomas_2)

# Create the globe bounding box
globe_bbox <- rbind(c(-180, -90), c(-180, 90), 
                    c(180, 90), c(180, -90), c(-180, -90))

# Create the globe border
globe_border_CAPTAIN2_EDGE2 <- st_polygon(list(globe_bbox)) %>%
  st_sfc(crs = 4326) %>%
  st_sf(data.frame(rgn = 'globe', geom = .)) %>%
  smoothr::densify(max_distance = 0.5) %>%
  st_transform(crs = mcbryde_thomas_2)

# Create base theme
my_theme_CAPTAIN2_EDGE2 <- theme_minimal() +
  theme(
    legend.position = "bottom",
    legend.direction = "horizontal",
    legend.box = "vertical",
    legend.margin = margin(t = 20, r = 0, b = 0, l = 0),
    legend.title = element_text(margin = margin(b = 10)),
    panel.background = element_rect(fill = "white", color = NA),
    plot.background = element_rect(fill = "white", color = NA),
    panel.grid = element_blank(),
    axis.text = element_blank(),      # Add this line
    axis.ticks = element_blank()      # Add this line
  )

# Create the plot
CAPTAIN2_EDGE2_plot <- ggplot() +
  geom_sf(data = CAPTAIN2_EDGE2_sf, aes(color = Priority), size = 0.5, alpha = 0.7) +
  geom_sf(data = world_projected_CAPTAIN2_EDGE2, fill = "lightgrey", color = "lightgrey", size = 0.1) +
  geom_sf(data = globe_border_CAPTAIN2_EDGE2, fill = NA, color = "black", size = 0.5) +
  scale_color_gradientn(
    colors = c("white", "yellow", "darkblue"),
    values = c(0, 0.5, 1),
    name = "Priority",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(title = "Global Conservation Priorities",
       subtitle = "CAPTAIN2 - EDGE2 Index, Budget: 0.1, Replicates: 50",
       x = NULL, y = NULL) +
  my_theme_CAPTAIN2_EDGE2

# Display the plot
print(CAPTAIN2_EDGE2_plot)

# Save the plot
ggsave(
  filename = here::here("Outputs", "CAPTAIN2_2ndrun_EDGE2_priorities_01_2.png"),
  plot = CAPTAIN2_EDGE2_plot,
  width = 10,
  height = 6,
  dpi = 300,
  bg = "white"
)

FUSE continental 0.1 budget

# Get the dimensions of the raster
nrows <- nrow(r)
ncols <- ncol(r)

# Confirm dimensions match expected values
if (nrows != 323 || ncols != 720) {
  warning("Raster dimensions don't match expected values. Proceeding with actual dimensions.")
}

# Create a grid of coordinates for each cell
# This gives us the center coordinates of each cell
coords <- as.data.frame(coordinates(r))
names(coords) <- c("Longitude", "Latitude")

# Add cell IDs (PUID) to the coordinates
coords$PUID <- 1:nrow(coords)

# Now join with the CAPTAIN2 FUSE data based on PUID
CAPTAIN2_FUSE_data_with_coords <- CAPTAIN2_FUSE_data %>%
  left_join(coords, by = "PUID")

# Check if the join worked correctly
if (sum(is.na(CAPTAIN2_FUSE_data_with_coords$Longitude)) > 0) {
  warning("Some PUIDs from CAPTAIN2 FUSE data couldn't be matched to coordinates.")
}

# Filter to keep only cells with non-zero priority for faster plotting
CAPTAIN2_FUSE_data_nonzero <- CAPTAIN2_FUSE_data_with_coords %>%
  filter(Priority > 0) %>%
  filter(!is.na(Longitude), !is.na(Latitude))  # Remove any rows with missing coords

# Get world map data
world <- ne_countries(scale = "medium", returnclass = "sf")

# Define the McBryde-Thomas 2 projection
mcbryde_thomas_2 <- "+proj=mbt_s"

# Transform the dataset to sf object and project
CAPTAIN2_FUSE_sf <- st_as_sf(
  CAPTAIN2_FUSE_data_nonzero, 
  coords = c("Longitude", "Latitude"), 
  crs = crs(r, asText = TRUE)  # Use the raster's CRS
) %>%
  st_transform(crs = mcbryde_thomas_2)

# Project the world map
world_projected_CAPTAIN2_FUSE <- st_transform(world, crs = mcbryde_thomas_2)

# Create the globe bounding box
globe_bbox <- rbind(c(-180, -90), c(-180, 90), 
                    c(180, 90), c(180, -90), c(-180, -90))

# Create the globe border
globe_border_CAPTAIN2_FUSE <- st_polygon(list(globe_bbox)) %>%
  st_sfc(crs = 4326) %>%
  st_sf(data.frame(rgn = 'globe', geom = .)) %>%
  smoothr::densify(max_distance = 0.5) %>%
  st_transform(crs = mcbryde_thomas_2)

# Create base theme
my_theme_CAPTAIN2_FUSE <- theme_minimal() +
  theme(
    legend.position = "bottom",
    legend.direction = "horizontal",
    legend.box = "vertical",
    legend.margin = margin(t = 20, r = 0, b = 0, l = 0),
    legend.title = element_text(margin = margin(b = 10)),
    panel.background = element_rect(fill = "white", color = NA),
    plot.background = element_rect(fill = "white", color = NA),
    panel.grid = element_blank(),
    axis.text = element_blank(),      # Add this line
    axis.ticks = element_blank())      # Add this line

# Create the plot
CAPTAIN2_FUSE_plot <- ggplot() +
  geom_sf(data = CAPTAIN2_FUSE_sf, aes(color = Priority), size = 0.5, alpha = 0.7) +
  geom_sf(data = world_projected_CAPTAIN2_FUSE, fill = "lightgrey", color = "lightgrey", size = 0.1) +
  geom_sf(data = globe_border_CAPTAIN2_FUSE, fill = NA, color = "black", size = 0.5) +
  scale_color_gradientn(
    colors = c("white", "yellow", "darkblue"),
    values = c(0, 0.5, 1),
    name = "Priority",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(title = "Global Conservation Priorities",
       subtitle = "CAPTAIN2 - FUSE Index, Budget: 0.1, Replicates: 50",
       x = NULL, y = NULL) +
  my_theme_CAPTAIN2_FUSE

# Display the plot
print(CAPTAIN2_FUSE_plot)

# Save the plot
ggsave(
  filename = here::here("Outputs", "CAPTAIN2_FUSE_priorities_01_2.png"),
  plot = CAPTAIN2_FUSE_plot,
  width = 10,
  height = 6,
  dpi = 300,
  bg = "white"
)

IUCN continental 0.1 budget

# Get the dimensions of the raster
nrows <- nrow(r)
ncols <- ncol(r)

# Confirm dimensions match expected values
if (nrows != 323 || ncols != 720) {
  warning("Raster dimensions don't match expected values. Proceeding with actual dimensions.")
}

# Create a grid of coordinates for each cell
# This gives us the center coordinates of each cell
coords <- as.data.frame(coordinates(r))
names(coords) <- c("Longitude", "Latitude")

# Add cell IDs (PUID) to the coordinates
coords$PUID <- 1:nrow(coords)

# Now join with the CAPTAIN2 IUCN data based on PUID
CAPTAIN2_IUCN_data_with_coords <- CAPTAIN2_IUCN_data %>%
  left_join(coords, by = "PUID")

# Check if the join worked correctly
if (sum(is.na(CAPTAIN2_IUCN_data_with_coords$Longitude)) > 0) {
  warning("Some PUIDs from CAPTAIN2 IUCN data couldn't be matched to coordinates.")
}

# Filter to keep only cells with non-zero priority for faster plotting
CAPTAIN2_IUCN_data_nonzero <- CAPTAIN2_IUCN_data_with_coords %>%
  filter(Priority > 0) %>%
  filter(!is.na(Longitude), !is.na(Latitude))  # Remove any rows with missing coords

# Get world map data
world <- ne_countries(scale = "medium", returnclass = "sf")

# Define the McBryde-Thomas 2 projection
mcbryde_thomas_2 <- "+proj=mbt_s"

# Transform the dataset to sf object and project
CAPTAIN2_IUCN_sf <- st_as_sf(
  CAPTAIN2_IUCN_data_nonzero, 
  coords = c("Longitude", "Latitude"), 
  crs = crs(r, asText = TRUE)  # Use the raster's CRS
) %>%
  st_transform(crs = mcbryde_thomas_2)

# Project the world map
world_projected_CAPTAIN2_IUCN <- st_transform(world, crs = mcbryde_thomas_2)

# Create the globe bounding box
globe_bbox <- rbind(c(-180, -90), c(-180, 90), 
                    c(180, 90), c(180, -90), c(-180, -90))

# Create the globe border
globe_border_CAPTAIN2_IUCN <- st_polygon(list(globe_bbox)) %>%
  st_sfc(crs = 4326) %>%
  st_sf(data.frame(rgn = 'globe', geom = .)) %>%
  smoothr::densify(max_distance = 0.5) %>%
  st_transform(crs = mcbryde_thomas_2)

# Create base theme
my_theme_CAPTAIN2_IUCN <- theme_minimal() +
  theme(
    legend.position = "bottom",
    legend.direction = "horizontal",
    legend.box = "vertical",
    legend.margin = margin(t = 20, r = 0, b = 0, l = 0),
    legend.title = element_text(margin = margin(b = 10)),
    panel.background = element_rect(fill = "white", color = NA),
    plot.background = element_rect(fill = "white", color = NA),
    panel.grid = element_blank(),
    axis.text = element_blank(),      # Add this line
    axis.ticks = element_blank()      # Add this line
  )

# Create the plot
CAPTAIN2_IUCN_plot <- ggplot() +
  geom_sf(data = CAPTAIN2_IUCN_sf, aes(color = Priority), size = 0.5, alpha = 0.7) +
  geom_sf(data = world_projected_CAPTAIN2_IUCN, fill = "lightgrey", color = "lightgrey", size = 0.1) +
  geom_sf(data = globe_border_CAPTAIN2_IUCN, fill = NA, color = "black", size = 0.5) +
  scale_color_gradientn(
    colors = c("white", "yellow", "darkblue"),
    values = c(0, 0.5, 1),
    name = "Priority",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(title = "Global Conservation Priorities",
       subtitle = "CAPTAIN2 - IUCN Index, Budget: 0.1, Replicates: 50",
       x = NULL, y = NULL) +
  my_theme_CAPTAIN2_IUCN

# Display the plot
print(CAPTAIN2_IUCN_plot)

# Save the plot
ggsave(
  filename = here::here("Outputs", "CAPTAIN2_IUCN_priorities_01_2.png"),
  plot = CAPTAIN2_IUCN_plot,
  width = 10,
  height = 6,
  dpi = 300,
  bg = "white"
)

Difference maps

# Create a grid of coordinates for each cell
coords <- as.data.frame(coordinates(r))
names(coords) <- c("Longitude", "Latitude")

# Add cell IDs (PUID) to the coordinates
coords$PUID <- 1:nrow(coords)

# Join all three datasets with coordinates
IUCN_with_coords <- CAPTAIN2_IUCN_data %>%
  dplyr::select(PUID, IUCN = Priority) %>%
  left_join(coords, by = "PUID")

EDGE2_with_coords <- CAPTAIN2_EDGE2_data %>%
  dplyr::select(PUID, EDGE2 = Priority) %>%
  left_join(coords, by = "PUID") 

FUSE_with_coords <- CAPTAIN2_FUSE_data %>%
  dplyr::select(PUID, FUSE = Priority) %>%
  left_join(coords, by = "PUID")

# Combine all datasets
all_indices <- coords %>%
  left_join(CAPTAIN2_IUCN_data %>% dplyr::select(PUID, IUCN = Priority), by = "PUID") %>%
  left_join(CAPTAIN2_EDGE2_data %>% dplyr::select(PUID, EDGE2 = Priority), by = "PUID") %>%
  left_join(CAPTAIN2_FUSE_data %>% dplyr::select(PUID, FUSE = Priority), by = "PUID")

# Calculate differences
all_indices <- all_indices %>%
  mutate(
    # Replace NA with 0 for calculation purposes
    IUCN = ifelse(is.na(IUCN), 0, IUCN),
    EDGE2 = ifelse(is.na(EDGE2), 0, EDGE2),
    FUSE = ifelse(is.na(FUSE), 0, FUSE),
    
    # Calculate differences
    IUCN_minus_FUSE = IUCN - FUSE,
    IUCN_minus_EDGE2 = IUCN - EDGE2,
    EDGE2_minus_FUSE = EDGE2 - FUSE
  )

# Get world map data
world <- ne_countries(scale = "medium", returnclass = "sf")

# Define the McBryde-Thomas 2 projection
mcbryde_thomas_2 <- "+proj=mbt_s"

# Create the globe bounding box
globe_bbox <- rbind(c(-180, -90), c(-180, 90), 
                    c(180, 90), c(180, -90), c(-180, -90))

# Create the globe border
globe_border <- st_polygon(list(globe_bbox)) %>%
  st_sfc(crs = 4326) %>%
  st_sf(data.frame(rgn = 'globe', geom = .)) %>%
  smoothr::densify(max_distance = 0.5) %>%
  st_transform(crs = mcbryde_thomas_2)

# Project the world map
world_projected <- st_transform(world, crs = mcbryde_thomas_2)

# Create base theme
my_theme <- theme_minimal() +
  theme(
    legend.position = "bottom",
    legend.direction = "horizontal",
    legend.box = "vertical",
    legend.margin = margin(t = 20, r = 0, b = 0, l = 0),
    legend.title = element_text(margin = margin(b = 10)),
    panel.background = element_rect(fill = "white", color = NA),
    plot.background = element_rect(fill = "white", color = NA),
    panel.grid = element_blank(),
    axis.text = element_blank(),      # Add this line
    axis.ticks = element_blank()      # Add this line
  )

# Filter to non-zero differences for each comparison to reduce plot size
# IUCN - FUSE
IUCN_FUSE_diff <- all_indices %>%
  filter(IUCN_minus_FUSE != 0) %>%
  st_as_sf(coords = c("Longitude", "Latitude"), crs = crs(r, asText = TRUE)) %>%
  st_transform(crs = mcbryde_thomas_2)

# IUCN - EDGE2
IUCN_EDGE2_diff <- all_indices %>%
  filter(IUCN_minus_EDGE2 != 0) %>%
  st_as_sf(coords = c("Longitude", "Latitude"), crs = crs(r, asText = TRUE)) %>%
  st_transform(crs = mcbryde_thomas_2)

# EDGE2 - FUSE
EDGE2_FUSE_diff <- all_indices %>%
  filter(EDGE2_minus_FUSE != 0) %>%
  st_as_sf(coords = c("Longitude", "Latitude"), crs = crs(r, asText = TRUE)) %>%
  st_transform(crs = mcbryde_thomas_2)

# Create a diverging color palette for difference maps
# Blue for negative (first index lower), white for zero, red for positive (first index higher)
diff_colors <- c("#2166AC", "#4393C3", "#92C5DE", "#D1E5F0", "#FFFFFF", 
                "#FDDBC7", "#F4A582", "#D6604D", "#B2182B")

# 1. IUCN - FUSE Difference Map
IUCN_FUSE_plot <- ggplot() +
  geom_sf(data = globe_border, fill = "#F8F8F8", color = NA) +
  geom_sf(data = IUCN_FUSE_diff, aes(color = IUCN_minus_FUSE), size = 0.5) +
  geom_sf(data = world_projected, fill = "lightgrey", color = "darkgrey", size = 0.1) +
  geom_sf(data = globe_border, fill = NA, color = "black", size = 0.5) +
  scale_color_gradientn(
    colors = diff_colors,
    limits = c(-1, 1),
    breaks = seq(-1, 1, by = 0.25),
    labels = as.character(seq(-1, 1, by = 0.25)),
    name = "Difference in Priority (IUCN - FUSE)",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(title = "Difference in Conservation Priorities",
       subtitle = "IUCN Index minus FUSE Index",
       x = NULL, y = NULL) +
  my_theme

# 2. IUCN - EDGE2 Difference Map
IUCN_EDGE2_plot <- ggplot() +
  geom_sf(data = globe_border, fill = "#F8F8F8", color = NA) +
  geom_sf(data = IUCN_EDGE2_diff, aes(color = IUCN_minus_EDGE2), size = 0.5) +
  geom_sf(data = world_projected, fill = "lightgrey", color = "darkgrey", size = 0.1) +
  geom_sf(data = globe_border, fill = NA, color = "black", size = 0.5) +
  scale_color_gradientn(
    colors = diff_colors,
    limits = c(-1, 1),
    breaks = seq(-1, 1, by = 0.25),
    labels = as.character(seq(-1, 1, by = 0.25)),
    name = "Difference in Priority (IUCN - EDGE2)",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(title = "Difference in Conservation Priorities",
       subtitle = "IUCN Index minus EDGE2 Index",
       x = NULL, y = NULL) +
  my_theme

# 3. EDGE2 - FUSE Difference Map
EDGE2_FUSE_plot <- ggplot() +
  geom_sf(data = globe_border, fill = "#F8F8F8", color = NA) +
  geom_sf(data = EDGE2_FUSE_diff, aes(color = EDGE2_minus_FUSE), size = 0.5) +
  geom_sf(data = world_projected, fill = "lightgrey", color = "darkgrey", size = 0.1) +
  geom_sf(data = globe_border, fill = NA, color = "black", size = 0.5) +
  scale_color_gradientn(
    colors = diff_colors,
    limits = c(-1, 1),
    breaks = seq(-1, 1, by = 0.25),
    labels = as.character(seq(-1, 1, by = 0.25)),
    name = "Difference in Priority (EDGE2 - FUSE)",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(title = "Difference in Conservation Priorities",
       subtitle = "EDGE2 Index minus FUSE Index",
       x = NULL, y = NULL) +
  my_theme

# Display all plots
print(IUCN_FUSE_plot)

print(IUCN_EDGE2_plot)

print(EDGE2_FUSE_plot)

# Save all plots
ggsave(
  filename = here::here("outputs", "CAPTAIN2_IUCN_minus_FUSE_difference_2.png"),
  plot = IUCN_FUSE_plot,
  width = 10,
  height = 6,
  dpi = 300,
  bg = "white"
)

ggsave(
  filename = here::here("outputs", "CAPTAIN2_IUCN_minus_EDGE2_difference_2.png"),
  plot = IUCN_EDGE2_plot,
  width = 10,
  height = 6,
  dpi = 300,
  bg = "white"
)

ggsave(
  filename = here::here("outputs", "CAPTAIN2_EDGE2_minus_FUSE_difference_2.png"),
  plot = EDGE2_FUSE_plot,
  width = 10,
  height = 6,
  dpi = 300,
  bg = "white"
)

Species level priorities

# Rename column in shark_metrics to match better
shark_metrics <- shark_metrics %>%
  rename(Species = `Species name`)

# Order shark_metrics alphabetically by Species name
shark_metrics <- shark_metrics %>%
  arrange(Species)

# Add an original order ID to protected_fractions to maintain its original order
protected_fractions$original_order <- 1:nrow(protected_fractions)

# Add row number as IDs to both datasets
protected_fractions$protected_ID <- 1:nrow(protected_fractions)
shark_metrics$species_ID <- 1:nrow(shark_metrics)

# Check if the datasets have the same number of rows
if(nrow(protected_fractions) == nrow(shark_metrics)) {
  # Create index mapping
  indices <- data.frame(
    protected_ID = 1:nrow(protected_fractions),
    species_ID = 1:nrow(shark_metrics)
  )
  
  # Join protected_fractions with indices
  protected_with_indices <- protected_fractions %>%
    left_join(indices, by = "protected_ID")
  
  # Join shark_metrics with indices
  shark_with_indices <- shark_metrics %>%
    left_join(indices, by = "species_ID")
  
  # Now join the datasets
  combined_data <- protected_with_indices %>%
    inner_join(
      shark_with_indices,
      by = c("species_ID", "protected_ID"),
      suffix = c("_captain", "_original")
    ) %>%
    arrange(original_order)
  
  cat("Successfully joined datasets with", nrow(combined_data), "species\n")
  
  # Define IUCN categories and order
  iucn_order <- c("LC", "NT", "VU", "EN", "CR")
  
  iucn_colors <- c(
    "LC" = "#50C878",
    "NT" = "#FFFF00",
    "VU" = "#FFA500",
    "EN" = "#FF8C00",
    "CR" = "#FF0000"
  )
  
  # Calculate sample sizes for each group
  iucn_n <- combined_data %>%
    group_by(IUCN_original) %>%
    summarise(n = n()) %>%
    arrange(IUCN_original)
  
  fuse_n <- combined_data %>%
    group_by(FUSE_original) %>%
    summarise(n = n()) %>%
    arrange(FUSE_original)
  
  edge2_n <- combined_data %>%
    group_by(EDGE2_original) %>%
    summarise(n = n()) %>%
    arrange(EDGE2_original)
  
  # Create x-axis labels with sample sizes
  iucn_labels <- paste0(iucn_order, "\n(n=", iucn_n$n, ")")
  fuse_labels <- paste0(1:5, "\n(n=", fuse_n$n, ")")
  edge2_labels <- paste0(1:5, "\n(n=", edge2_n$n, ")")
  
  # 1. IUCN plot
  iucn_plot <- combined_data %>%
    mutate(
      IUCN_status = factor(IUCN_original, levels = 1:5, labels = iucn_order),
      protection_percentage = IUCN_captain * 100
    ) %>%
    ggplot(aes(x = IUCN_status, y = protection_percentage)) +
    geom_jitter(width = 0.1,
                size = 0.4,
                alpha = 0.5,
                color = "darkgray") +
    stat_summary(fun = mean,
                 geom = "point",
                 size = 3,
                 color = "black") +
    stat_summary(fun.data = function(x) {
      mean_val <- mean(x, na.rm = TRUE)
      sd_val <- sd(x, na.rm = TRUE)
      return(data.frame(y = mean_val,
                        ymin = max(0, mean_val - sd_val),
                        ymax = min(100, mean_val + sd_val)))
    },
    geom = "errorbar",
    width = 0.1,
    color = "black") +
    labs(
      x = "IUCN Red List threat status",
      y = "Range protected (%)"
    ) +
    theme_minimal() +
    theme(
      plot.title = element_text(hjust = 0.5, face = "bold", size = 18),
      legend.position = "none",
      panel.grid.major.x = element_blank(),
      plot.background = element_rect(fill = "white", color = NA),
      panel.background = element_rect(fill = "white", color = NA),
      axis.title.x = element_text(size = 16),
      axis.title.y = element_text(size = 16),
      axis.text.x = element_text(size = 14),
      axis.text.y = element_text(size = 14)
    ) +
    scale_x_discrete(limits = iucn_order, labels = iucn_labels) +
    scale_fill_manual(values = iucn_colors) +
    scale_color_manual(values = iucn_colors) +
    scale_y_continuous(limits = c(0, 100),
                       breaks = seq(0, 100, 25))
  
  # 2. FUSE plot
  fuse_plot <- combined_data %>%
    mutate(
      FUSE_category = factor(FUSE_original),
      protection_percentage = FUSE_captain * 100
    ) %>%
    ggplot(aes(x = FUSE_category, y = protection_percentage)) +
    geom_jitter(width = 0.1,
                size = 0.4,
                alpha = 0.5,
                color = "darkgray") +
    stat_summary(fun = mean,
                 geom = "point",
                 size = 3,
                 color = "black") +
    stat_summary(fun.data = function(x) {
      mean_val <- mean(x, na.rm = TRUE)
      sd_val <- sd(x, na.rm = TRUE)
      return(data.frame(y = mean_val,
                        ymin = max(0, mean_val - sd_val),
                        ymax = min(100, mean_val + sd_val)))
    },
    geom = "errorbar",
    width = 0.1,
    color = "black") +
    labs(
      x = "FUSE Score",
      y = "Range protected (%)"
    ) +
    theme_minimal() +
    theme(
      plot.title = element_text(hjust = 0.5, face = "bold", size = 18),
      legend.position = "none",
      panel.grid.major.x = element_blank(),
      plot.background = element_rect(fill = "white", color = NA),
      panel.background = element_rect(fill = "white", color = NA),
      axis.title.x = element_text(size = 16),
      axis.title.y = element_text(size = 16),
      axis.text.x = element_text(size = 14),
      axis.text.y = element_text(size = 14)
    ) +
    scale_x_discrete(labels = fuse_labels) +
    scale_y_continuous(limits = c(0, 100),
                       breaks = seq(0, 100, 25))
  
  # 3. EDGE2 plot
  edge2_plot <- combined_data %>%
    mutate(
      EDGE2_category = factor(EDGE2_original),
      protection_percentage = EDGE2_captain * 100
    ) %>%
    ggplot(aes(x = EDGE2_category, y = protection_percentage)) +
    geom_jitter(width = 0.1,
                size = 0.4,
                alpha = 0.5,
                color = "darkgray") +
    stat_summary(fun = mean,
                 geom = "point",
                 size = 3,
                 color = "black") +
    stat_summary(fun.data = function(x) {
      mean_val <- mean(x, na.rm = TRUE)
      sd_val <- sd(x, na.rm = TRUE)
      return(data.frame(y = mean_val,
                        ymin = max(0, mean_val - sd_val),
                        ymax = min(100, mean_val + sd_val)))
    },
    geom = "errorbar",
    width = 0.1,
    color = "black") +
    labs(
      x = "EDGE2 Score",
      y = "Range protected (%)"
    ) +
    theme_minimal() +
    theme(
      plot.title = element_text(hjust = 0.5, face = "bold", size = 18),
      legend.position = "none",
      panel.grid.major.x = element_blank(),
      plot.background = element_rect(fill = "white", color = NA),
      panel.background = element_rect(fill = "white", color = NA),
      axis.title.x = element_text(size = 16),
      axis.title.y = element_text(size = 16),
      axis.text.x = element_text(size = 14),
      axis.text.y = element_text(size = 14)
    ) +
    scale_x_discrete(labels = edge2_labels) +
    scale_y_continuous(limits = c(0, 100),
                       breaks = seq(0, 100, 25))
  
  # Combine the three plots into a grid
  combined_plots <- ggpubr::ggarrange(
    iucn_plot,
    fuse_plot,
    edge2_plot,
    ncol = 3,
    nrow = 1,
    labels = c("(A)", "(B)", "(C)"),
    font.label = list(size = 14,
                      face = "bold",
                      hjust = -2,
                      vjust = 1.5)
  )
  
  print(combined_plots)
  
  # Save combined plot
  ggsave(here::here("Outputs", "Fig.2_combined_protection_dotplots.png"),
         combined_plots, width = 20, height = 6, dpi = 150, bg = "white")
  
  #======================
  # Beta Regression Analysis
  #======================
  
  # Prepare data for beta regression
  model_data <- combined_data %>%
    mutate(
      IUCN_status = IUCN_original,
      FUSE_score = FUSE_original,
      EDGE2_score = EDGE2_original,
      IUCN_protection = IUCN_captain * 100,
      FUSE_protection = FUSE_captain * 100,
      EDGE2_protection = EDGE2_captain * 100,
      # Convert to (0,1) scale and handle boundaries
      n_obs = n(),
      IUCN_protection_beta = (IUCN_protection/100 * (n_obs - 1) + 0.5) / n_obs,
      FUSE_protection_beta = (FUSE_protection/100 * (n_obs - 1) + 0.5) / n_obs,
      EDGE2_protection_beta = (EDGE2_protection/100 * (n_obs - 1) + 0.5) / n_obs
    )
  
  # Fit beta regression models
  iucn_beta <- betareg(IUCN_protection_beta ~ IUCN_status, data = model_data)
  fuse_beta <- betareg(FUSE_protection_beta ~ FUSE_score, data = model_data)
  edge2_beta <- betareg(EDGE2_protection_beta ~ EDGE2_score, data = model_data)
  
  # View summaries
  cat("\n=== IUCN Beta Regression ===\n")
  print(summary(iucn_beta))
  
  cat("\n=== FUSE Beta Regression ===\n")
  print(summary(fuse_beta))
  
  cat("\n=== EDGE2 Beta Regression ===\n")
  print(summary(edge2_beta))
  
  # Get tidy summaries
  iucn_tidy <- tidy(iucn_beta)
  fuse_tidy <- tidy(fuse_beta)
  edge2_tidy <- tidy(edge2_beta)
  
  # Get pseudo R²
  iucn_r2 <- cor(predict(iucn_beta), model_data$IUCN_protection_beta)^2
  fuse_r2 <- cor(predict(fuse_beta), model_data$FUSE_protection_beta)^2
  edge2_r2 <- cor(predict(edge2_beta), model_data$EDGE2_protection_beta)^2
  
  # Calculate marginal effects (average change in percentage per unit increase)
  iucn_margins <- margins(iucn_beta)
  fuse_margins <- margins(fuse_beta)
  edge2_margins <- margins(edge2_beta)
  
  # Extract AME and SE
  iucn_ame_summary <- summary(iucn_margins)
  fuse_ame_summary <- summary(fuse_margins)
  edge2_ame_summary <- summary(edge2_margins)
  
  # Convert to percentage (margins gives proportions)
  iucn_ame <- iucn_ame_summary$AME * 100
  iucn_ame_se <- iucn_ame_summary$SE * 100
  
  fuse_ame <- fuse_ame_summary$AME * 100
  fuse_ame_se <- fuse_ame_summary$SE * 100
  
  edge2_ame <- edge2_ame_summary$AME * 100
  edge2_ame_se <- edge2_ame_summary$SE * 100
  
  cat("\n=== Average Marginal Effects (as percentages) ===\n")
  cat(paste("IUCN AME:", round(iucn_ame, 2), "%, SE:", round(iucn_ame_se, 2), "%\n"))
  cat(paste("FUSE AME:", round(fuse_ame, 2), "%, SE:", round(fuse_ame_se, 2), "%\n"))
  cat(paste("EDGE2 AME:", round(edge2_ame, 2), "%, SE:", round(edge2_ame_se, 2), "%\n"))
  
  # Calculate predicted values and differences
  # IUCN predictions
  iucn_pred <- predict(iucn_beta,
                       newdata = data.frame(IUCN_status = 1:5),
                       type = "response") * 100
  
  cat("\nPredicted protection by IUCN category:\n")
  print(data.frame(Category = 1:5,
                   IUCN_Category = c("LC", "NT", "VU", "EN", "CR"),
                   Predicted_Protection = round(iucn_pred, 2)))
  
  iucn_avg_increase <- mean(diff(iucn_pred))
  cat(paste("Average increase per IUCN category:", round(iucn_avg_increase, 2), "%\n"))
  
  # FUSE predictions
  fuse_pred <- predict(fuse_beta,
                       newdata = data.frame(FUSE_score = 1:5),
                       type = "response") * 100
  
  cat("\nPredicted protection by FUSE score:\n")
  print(data.frame(Score = 1:5,
                   Predicted_Protection = round(fuse_pred, 2)))
  
  fuse_avg_increase <- mean(diff(fuse_pred))
  cat(paste("Average increase per FUSE unit:", round(fuse_avg_increase, 2), "%\n"))
  
  # EDGE2 predictions
  edge2_pred <- predict(edge2_beta,
                        newdata = data.frame(EDGE2_score = 1:5),
                        type = "response") * 100
  
  cat("\nPredicted protection by EDGE2 score:\n")
  print(data.frame(Score = 1:5,
                   Predicted_Protection = round(edge2_pred, 2)))
  
  edge2_avg_increase <- mean(diff(edge2_pred))
  cat(paste("Average increase per EDGE2 unit:", round(edge2_avg_increase, 2), "%\n"))
  
  # Create comprehensive summary table with AME
  comprehensive_summary <- data.frame(
    Approach = c("IUCN", "FUSE", "EDGE2"),
    Coefficient_logit = c(
      iucn_tidy$estimate[2],
      fuse_tidy$estimate[2],
      edge2_tidy$estimate[2]
    ),
    SE = c(
      iucn_tidy$std.error[2],
      fuse_tidy$std.error[2],
      edge2_tidy$std.error[2]
    ),
    z_value = c(
      iucn_tidy$statistic[2],
      fuse_tidy$statistic[2],
      edge2_tidy$statistic[2]
    ),
    P_value = c(
      iucn_tidy$p.value[2],
      fuse_tidy$p.value[2],
      edge2_tidy$p.value[2]
    ),
    AME_percent = c(
      iucn_ame,
      fuse_ame,
      edge2_ame
    ),
    AME_SE = c(
      iucn_ame_se,
      fuse_ame_se,
      edge2_ame_se
    ),
    Pseudo_R2 = c(iucn_r2, fuse_r2, edge2_r2)
  )
  
  cat("\n=== Comprehensive Beta Regression Summary ===\n")
  print(comprehensive_summary)
  
  # Save results
  write.csv(comprehensive_summary,
            here::here("Outputs", "beta_regression_summary.csv"),
            row.names = FALSE)
  
} else {
  cat("ERROR: Datasets have different number of rows.\n")
  cat("Protected fractions:", nrow(protected_fractions), "rows\n")
  cat("Shark metrics:", nrow(shark_metrics), "rows\n")
}

Successfully joined datasets with 1000 species

=== IUCN Beta Regression ===

Call:
betareg(formula = IUCN_protection_beta ~ IUCN_status, data = model_data)

Quantile residuals:
    Min      1Q  Median      3Q     Max 
-2.5313 -0.5004 -0.0020  0.5608  2.0202 

Coefficients (mean model with logit link):
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.05978    0.07511  -0.796    0.426    
IUCN_status  0.11768    0.02988   3.938 8.21e-05 ***

Phi coefficients (precision model with identity link):
      Estimate Std. Error z value Pr(>|z|)    
(phi)  0.91473    0.03196   28.63   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Type of estimator: ML (maximum likelihood)
Log-likelihood: 379.5 on 3 Df
Pseudo R-squared: 0.01504
Number of iterations: 7 (BFGS) + 1 (Fisher scoring) 

=== FUSE Beta Regression ===

Call:
betareg(formula = FUSE_protection_beta ~ FUSE_score, data = model_data)

Quantile residuals:
    Min      1Q  Median      3Q     Max 
-2.3375 -0.4546  0.0051  0.5057  2.3049 

Coefficients (mean model with logit link):
            Estimate Std. Error z value Pr(>|z|)  
(Intercept) -0.19664    0.08966  -2.193   0.0283 *
FUSE_score   0.14261    0.06830   2.088   0.0368 *

Phi coefficients (precision model with identity link):
      Estimate Std. Error z value Pr(>|z|)    
(phi)  1.04726    0.03673   28.52   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Type of estimator: ML (maximum likelihood)
Log-likelihood: 211.9 on 3 Df
Pseudo R-squared: 0.004321
Number of iterations: 7 (BFGS) + 2 (Fisher scoring) 

=== EDGE2 Beta Regression ===

Call:
betareg(formula = EDGE2_protection_beta ~ EDGE2_score, data = model_data)

Quantile residuals:
    Min      1Q  Median      3Q     Max 
-2.5490 -0.4954 -0.0052  0.4548  2.1662 

Coefficients (mean model with logit link):
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.68964    0.09615  -7.173 7.34e-13 ***
EDGE2_score  0.46406    0.07688   6.036 1.58e-09 ***

Phi coefficients (precision model with identity link):
      Estimate Std. Error z value Pr(>|z|)    
(phi)  0.84174    0.02926   28.77   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Type of estimator: ML (maximum likelihood)
Log-likelihood: 506.7 on 3 Df
Pseudo R-squared: 0.04761
Number of iterations: 13 (BFGS) + 1 (Fisher scoring) 

=== Average Marginal Effects (as percentages) ===
IUCN AME: 2.9 %, SE: 0.73 %
FUSE AME: 3.56 %, SE: 1.7 %
EDGE2 AME: 11.39 %, SE: 1.85 %

Predicted protection by IUCN category:
  Category IUCN_Category Predicted_Protection
1        1            LC                51.45
2        2            NT                54.38
3        3            VU                57.28
4        4            EN                60.13
5        5            CR                62.92
Average increase per IUCN category: 2.87 %

Predicted protection by FUSE score:
  Score Predicted_Protection
1     1                48.65
2     2                52.21
3     3                55.75
4     4                59.24
5     5                62.63
Average increase per FUSE unit: 3.5 %

Predicted protection by EDGE2 score:
  Score Predicted_Protection
1     1                44.38
2     2                55.93
3     3                66.88
4     4                76.25
5     5                83.63
Average increase per EDGE2 unit: 9.81 %

=== Comprehensive Beta Regression Summary ===
            Approach Coefficient_logit         SE  z_value      P_value
IUCN_status     IUCN         0.1176756 0.02987957 3.938332 8.205012e-05
FUSE_score      FUSE         0.1426108 0.06830262 2.087925 3.680458e-02
EDGE2_score    EDGE2         0.4640609 0.07687772 6.036350 1.576390e-09
            AME_percent    AME_SE  Pseudo_R2
IUCN_status    2.898230 0.7271499 0.01750165
FUSE_score     3.558461 1.6986270 0.00239858
EDGE2_score   11.386873 1.8514148 0.02982912

Manuscript maps

Individual maps

CAPTAIN2_EDGE2_msmap <- ggplot() +
  geom_sf(data = CAPTAIN2_EDGE2_sf, aes(color = Priority), size = 0.5, alpha = 0.7) +
  geom_sf(data = world_projected_CAPTAIN2_EDGE2, fill = "lightgrey", color = "lightgrey", size = 0.1) +
  geom_sf(data = globe_border_CAPTAIN2_EDGE2, fill = NA, color = "lightgrey", size = 0.5) +
  scale_color_gradientn(
    colors = c("white", "yellow", "darkblue"),
    values = c(0, 0.5, 1),
    name = "Priority EDGE2",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(#title = "Global Conservation Priorities",
       #subtitle = "CAPTAIN2 - EDGE2 Index, Budget: 0.1, Replicates: 50",
       x = NULL, y = NULL) +
  my_theme_CAPTAIN2_EDGE2
CAPTAIN2_EDGE2_msmap

CAPTAIN2_FUSE_msmap <- ggplot() +
  geom_sf(data = CAPTAIN2_FUSE_sf, aes(color = Priority), size = 0.5, alpha = 0.7) +
  geom_sf(data = world_projected_CAPTAIN2_FUSE, fill = "lightgrey", color = "lightgrey", size = 0.1) +
  geom_sf(data = globe_border_CAPTAIN2_FUSE, fill = NA, color = "lightgrey", size = 0.5) +
  scale_color_gradientn(
    colors = c("white", "yellow", "darkblue"),
    values = c(0, 0.5, 1),
    name = "Priority FUSE",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(#title = "Global Conservation Priorities",
       #subtitle = "CAPTAIN2 - FUSE Index, Budget: 0.1, Replicates: 50",
       x = NULL, y = NULL) +
  my_theme_CAPTAIN2_FUSE
CAPTAIN2_FUSE_msmap

CAPTAIN2_IUCN_msmap <- ggplot() +
  geom_sf(data = CAPTAIN2_IUCN_sf, aes(color = Priority), size = 0.5, alpha = 0.7) +
  geom_sf(data = world_projected_CAPTAIN2_IUCN, fill = "lightgrey", color = "lightgrey", size = 0.1) +
  geom_sf(data = globe_border_CAPTAIN2_IUCN, fill = NA, color = "lightgrey", size = 0.5) +
  scale_color_gradientn(
    colors = c("white", "yellow", "darkblue"),
    values = c(0, 0.5, 1),
    name = "Priority IUCN",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(#title = "Global Conservation Priorities",
       #subtitle = "CAPTAIN2 - IUCN Index, Budget: 0.1, Replicates: 50",
       x = NULL, y = NULL) +
  my_theme_CAPTAIN2_IUCN
CAPTAIN2_IUCN_msmap

# Create a function to add labels (A), (B), etc.
add_panel_labels <- function(plot, label) {
  plot + 
    theme(
      plot.title = element_text(face = "bold", hjust = 0, size = 12)
    ) +
    labs(title = paste0("(", label, ")"))
}

# Add labels to each plot
# First grid
CAPTAIN2_IUCN_msmap_labeled <- add_panel_labels(CAPTAIN2_IUCN_msmap, "A")
CAPTAIN2_FUSE_msmap_labeled <- add_panel_labels(CAPTAIN2_FUSE_msmap, "B")
CAPTAIN2_EDGE2_msmap_labeled <- add_panel_labels(CAPTAIN2_EDGE2_msmap, "C")

# Combine plots into two separate grids, each with 3 rows and 1 column
grid1 <- CAPTAIN2_IUCN_msmap_labeled /
         CAPTAIN2_FUSE_msmap_labeled /
         CAPTAIN2_EDGE2_msmap_labeled

# Save the grids if needed
ggsave(here::here("Outputs", "grid1_maps_continental_2ndrun_new.png"), grid1, width = 8, height = 15, dpi = 300)

Difference maps

# 1. IUCN - FUSE Difference Map
IUCN_FUSE_msmap <- ggplot() +
  geom_sf(data = globe_border, fill = "#F8F8F8", color = NA) +
  geom_sf(data = IUCN_FUSE_diff, aes(color = IUCN_minus_FUSE), size = 0.5) +
  geom_sf(data = world_projected, fill = "lightgrey", color = "darkgrey", size = 0.1) +
  geom_sf(data = globe_border, fill = NA, color = "lightgrey", size = 0.5) +
  scale_color_gradientn(
    colors = diff_colors,
    limits = c(-1, 1),
    breaks = seq(-1, 1, by = 0.25),
    labels = as.character(seq(-1, 1, by = 0.25)),
    name = "Difference in Priority (IUCN - FUSE)",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(#title = "Difference in Conservation Priorities",
       #subtitle = "IUCN Index minus FUSE Index",
       x = NULL, y = NULL) +
  my_theme

# 2. IUCN - EDGE2 Difference Map
IUCN_EDGE2_msmap <- ggplot() +
  geom_sf(data = globe_border, fill = "#F8F8F8", color = NA) +
  geom_sf(data = IUCN_EDGE2_diff, aes(color = IUCN_minus_EDGE2), size = 0.5) +
  geom_sf(data = world_projected, fill = "lightgrey", color = "darkgrey", size = 0.1) +
  geom_sf(data = globe_border, fill = NA, color = "lightgrey", size = 0.5) +
  scale_color_gradientn(
    colors = diff_colors,
    limits = c(-1, 1),
    breaks = seq(-1, 1, by = 0.25),
    labels = as.character(seq(-1, 1, by = 0.25)),
    name = "Difference in Priority (IUCN - EDGE2)",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(#title = "Difference in Conservation Priorities",
       #subtitle = "IUCN Index minus EDGE2 Index",
       x = NULL, y = NULL) +
  my_theme

# 3. EDGE2 - FUSE Difference Map
EDGE2_FUSE_msmap <- ggplot() +
  geom_sf(data = globe_border, fill = "#F8F8F8", color = NA) +
  geom_sf(data = EDGE2_FUSE_diff, aes(color = EDGE2_minus_FUSE), size = 0.5) +
  geom_sf(data = world_projected, fill = "lightgrey", color = "darkgrey", size = 0.1) +
  geom_sf(data = globe_border, fill = NA, color = "lightgrey", size = 0.5) +
  scale_color_gradientn(
    colors = diff_colors,
    limits = c(-1, 1),
    breaks = seq(-1, 1, by = 0.25),
    labels = as.character(seq(-1, 1, by = 0.25)),
    name = "Difference in Priority (EDGE2 - FUSE)",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(#title = "Difference in Conservation Priorities",
       #subtitle = "EDGE2 Index minus FUSE Index",
       x = NULL, y = NULL) +
  my_theme

# Create a function to add labels (A), (B), etc.
add_panel_labels <- function(plot, label) {
  plot + 
    theme(
      plot.title = element_text(face = "bold", hjust = 0, size = 12)
    ) +
    labs(title = paste0("(", label, ")"))
}

# Add labels to each plot
# First grid
IUCN_FUSE_msmap_labeled <- add_panel_labels(IUCN_FUSE_msmap, "A")
IUCN_EDGE2_msmap_labeled <- add_panel_labels(IUCN_EDGE2_msmap, "B")
EDGE2_FUSE_msmap_labeled <- add_panel_labels(EDGE2_FUSE_msmap, "C")

# Combine plots into two separate grids, each with 3 rows and 1 column
grid2 <- IUCN_FUSE_msmap_labeled /
         IUCN_EDGE2_msmap_labeled /
         EDGE2_FUSE_msmap_labeled

# Save the grids if needed
ggsave(here::here("Outputs", "grid3_maps_continental_2ndrun.png"), grid2, width = 16, height = 15, dpi = 300)

Combined individual and difference maps

library(patchwork)

# Add labels to each plot
# First grid
CAPTAIN2_IUCN_msmap_labeled <- add_panel_labels(CAPTAIN2_IUCN_msmap, "A")
CAPTAIN2_FUSE_msmap_labeled <- add_panel_labels(CAPTAIN2_FUSE_msmap, "B")
CAPTAIN2_EDGE2_msmap_labeled <- add_panel_labels(CAPTAIN2_EDGE2_msmap, "C")

# Second grid
IUCN_FUSE_msmap_labeled <- add_panel_labels(IUCN_FUSE_msmap, "D")
IUCN_EDGE2_msmap_labeled <- add_panel_labels(IUCN_EDGE2_msmap, "E")
EDGE2_FUSE_msmap_labeled <- add_panel_labels(EDGE2_FUSE_msmap, "F")

# Create each grid (3 rows, 1 column)
grid1 <- CAPTAIN2_IUCN_msmap_labeled /
         CAPTAIN2_FUSE_msmap_labeled /
         CAPTAIN2_EDGE2_msmap_labeled

grid2 <- IUCN_FUSE_msmap_labeled /
         IUCN_EDGE2_msmap_labeled /
         EDGE2_FUSE_msmap_labeled

# Combine the two grids side by side (2 columns)
combined_grid <- grid1 | grid2

# Display the combined grid
combined_grid

# If you want to save it
ggsave(
  filename = here::here("Outputs", "Fig.3_combined_priority_maps.png"),
  plot = combined_grid,
  width = 10,    # Adjust width as needed for two columns
  height = 12,   # Adjust height as needed for three rows
  dpi = 300,
  bg = "white"
)

Combined individual and difference maps

library(patchwork)

# Add labels to each plot
# First grid
CAPTAIN2_IUCN_msmap_labeled <- add_panel_labels(CAPTAIN2_IUCN_msmap, "A")
CAPTAIN2_FUSE_msmap_labeled <- add_panel_labels(CAPTAIN2_FUSE_msmap, "B")
CAPTAIN2_EDGE2_msmap_labeled <- add_panel_labels(CAPTAIN2_EDGE2_msmap, "C")

# Second grid
IUCN_FUSE_msmap_labeled <- add_panel_labels(IUCN_FUSE_msmap, "D")
IUCN_EDGE2_msmap_labeled <- add_panel_labels(IUCN_EDGE2_msmap, "E")
EDGE2_FUSE_msmap_labeled <- add_panel_labels(EDGE2_FUSE_msmap, "F")

# Create each grid (3 rows, 1 column)
grid1 <- CAPTAIN2_IUCN_msmap_labeled /
         CAPTAIN2_FUSE_msmap_labeled /
         CAPTAIN2_EDGE2_msmap_labeled

grid2 <- IUCN_FUSE_msmap_labeled /
         IUCN_EDGE2_msmap_labeled /
         EDGE2_FUSE_msmap_labeled

# Combine the two grids side by side (2 columns)
combined_grid <- grid1 | grid2

# Display the combined grid
combined_grid

# If you want to save it
ggsave(
  filename = here::here("outputs", "combined_priority_maps.png"),
  plot = combined_grid,
  width = 10,    # Adjust width as needed for two columns
  height = 12,   # Adjust height as needed for three rows
  dpi = 300,
  bg = "white"
)

Spatial distribution analysis

# Filter for high priority cells (>0.9) and calculate percentages
# Calculate total coastal cells and high-priority cells for each index
# IUCN
CAPTAIN2_IUCN_high <- CAPTAIN2_IUCN_sf %>% 
  filter(Priority > 0.9)
iucn_total_cells <- nrow(CAPTAIN2_IUCN_sf)
iucn_high_cells <- nrow(CAPTAIN2_IUCN_high)
iucn_percent <- round((iucn_high_cells / iucn_total_cells) * 100, 1)

# FUSE
CAPTAIN2_FUSE_high <- CAPTAIN2_FUSE_sf %>% 
  filter(Priority > 0.9)
fuse_total_cells <- nrow(CAPTAIN2_FUSE_sf)
fuse_high_cells <- nrow(CAPTAIN2_FUSE_high)
fuse_percent <- round((fuse_high_cells / fuse_total_cells) * 100, 1)

# EDGE2
CAPTAIN2_EDGE2_high <- CAPTAIN2_EDGE2_sf %>% 
  filter(Priority > 0.9)
edge2_total_cells <- nrow(CAPTAIN2_EDGE2_sf)
edge2_high_cells <- nrow(CAPTAIN2_EDGE2_high)
edge2_percent <- round((edge2_high_cells / edge2_total_cells) * 100, 1)

# Print results
cat("IUCN: ", iucn_high_cells, " cells (", iucn_percent, "%) with priority >0.9\n", sep="")

IUCN: 4734 cells (78.5%) with priority >0.9

cat("FUSE: ", fuse_high_cells, " cells (", fuse_percent, "%) with priority >0.9\n", sep="")

FUSE: 4396 cells (60.5%) with priority >0.9

cat("EDGE2: ", edge2_high_cells, " cells (", edge2_percent, "%) with priority >0.9\n", sep="")

EDGE2: 5328 cells (94.8%) with priority >0.9

#Mean nearest neighbor for high priority grid cells:

# Function to calculate Mean Nearest Neighbor Distance
calculate_mnn <- function(sf_high_priority) {
  # Get centroids and convert to coordinates
  coords <- st_coordinates(st_centroid(sf_high_priority))
  
  # Create a bounding box for the window
  bbox <- st_bbox(sf_high_priority)
  
  # Create observation window
  window <- owin(xrange = c(bbox["xmin"], bbox["xmax"]),
                 yrange = c(bbox["ymin"], bbox["ymax"]))
  
  # Create point pattern object
  ppp_obj <- ppp(coords[,1], coords[,2], window = window)
  
  # Calculate mean nearest neighbor distance
  mnn <- mean(nndist(ppp_obj))
  
  # Calculate standard deviation for context
  mnn_sd <- sd(nndist(ppp_obj))
  
  return(list(mean = mnn, sd = mnn_sd))
}

# Calculate for each index
iucn_mnn <- calculate_mnn(CAPTAIN2_IUCN_high)
fuse_mnn <- calculate_mnn(CAPTAIN2_FUSE_high)
edge2_mnn <- calculate_mnn(CAPTAIN2_EDGE2_high)

# Print results
cat("\nMean Nearest Neighbor Distance (units depend on your projection):\n")

Mean Nearest Neighbor Distance (units depend on your projection):

cat("IUCN:  Mean = ", round(iucn_mnn$mean, 2), ", SD = ", round(iucn_mnn$sd, 2), "\n", sep="")

IUCN:  Mean = 51143.37, SD = 17182.87

cat("FUSE:  Mean = ", round(fuse_mnn$mean, 2), ", SD = ", round(fuse_mnn$sd, 2), "\n", sep="")

FUSE:  Mean = 49594.94, SD = 26765.94

cat("EDGE2: Mean = ", round(edge2_mnn$mean, 2), ", SD = ", round(edge2_mnn$sd, 2), "\n", sep="")

EDGE2: Mean = 48211.63, SD = 8975.76

# Create summary data frame
mnn_summary <- data.frame(
  Index = c("IUCN", "FUSE", "EDGE2"),
  N_cells = c(iucn_high_cells, fuse_high_cells, edge2_high_cells),
  Percent = c(iucn_percent, fuse_percent, edge2_percent),
  MNN_mean = c(iucn_mnn$mean, fuse_mnn$mean, edge2_mnn$mean),
  MNN_sd = c(iucn_mnn$sd, fuse_mnn$sd, edge2_mnn$sd)
)

print(mnn_summary)

  Index N_cells Percent MNN_mean    MNN_sd
1  IUCN    4734    78.5 51143.37 17182.874
2  FUSE    4396    60.5 49594.94 26765.936
3 EDGE2    5328    94.8 48211.63  8975.762

Congruence maps

# Method 1: If you have separate dataframes for each index
# Assuming your data has columns: PUID, Priority, and geometry

# Fxirst, identify high priority areas (>0.9) for each index
high_priority_EDGE2 <- CAPTAIN2_EDGE2_sf[CAPTAIN2_EDGE2_sf$Priority > 0.9, ]
high_priority_FUSE <- CAPTAIN2_FUSE_sf[CAPTAIN2_FUSE_sf$Priority > 0.9, ]
high_priority_IUCN <- CAPTAIN2_IUCN_sf[CAPTAIN2_IUCN_sf$Priority > 0.9, ]

# Find congruent areas (present in all three)
# Method using PUID (assuming you have PUID column)
congruent_PUIDs <- intersect(intersect(high_priority_EDGE2$PUID, 
                                      high_priority_FUSE$PUID), 
                            high_priority_IUCN$PUID)

# Extract congruent areas
congruent_areas <- CAPTAIN2_EDGE2_sf[CAPTAIN2_EDGE2_sf$PUID %in% congruent_PUIDs, ]

# Print summary
cat("High priority areas (>0.9):\n")

High priority areas (>0.9):

cat("EDGE2:", nrow(high_priority_EDGE2), "areas\n")

EDGE2: 5328 areas

cat("FUSE:", nrow(high_priority_FUSE), "areas\n") 

FUSE: 4396 areas

cat("IUCN:", nrow(high_priority_IUCN), "areas\n")

IUCN: 4734 areas

cat("Congruent areas:", nrow(congruent_areas), "areas\n")

Congruent areas: 1459 areas

cat("Percentage of overlap:", round(nrow(congruent_areas)/min(nrow(high_priority_EDGE2), 
                                                            nrow(high_priority_FUSE), 
                                                            nrow(high_priority_IUCN))*100, 2), "%\n")

Percentage of overlap: 33.19 %

# Method 2: If you need to merge dataframes first
# Create a combined dataframe with all three priority scores
# First, extract just the data (without geometry) from the other sf objects
FUSE_data <- st_drop_geometry(CAPTAIN2_FUSE_sf[, c("PUID", "Priority")])
IUCN_data <- st_drop_geometry(CAPTAIN2_IUCN_sf[, c("PUID", "Priority")])

# Merge with the EDGE2 sf object (keeping geometry)
combined_priorities <- merge(CAPTAIN2_EDGE2_sf, FUSE_data, 
                            by = "PUID", suffixes = c("_EDGE2", "_FUSE"))
combined_priorities <- merge(combined_priorities, IUCN_data, 
                            by = "PUID")
names(combined_priorities)[names(combined_priorities) == "Priority"] <- "Priority_IUCN"

# Identify congruent high priority areas
congruent_areas_v2 <- combined_priorities[combined_priorities$Priority_EDGE2 > 0.9 & 
                                         combined_priorities$Priority_FUSE > 0.9 & 
                                         combined_priorities$Priority_IUCN > 0.9, ]

# Create a map showing only congruent areas
congruent_map <- ggplot() +
  geom_sf(data = congruent_areas, aes(fill = "Congruent High Priority"), 
          color = "red", size = 0.8, alpha = 0.8) +
  geom_sf(data = world_projected_CAPTAIN2_EDGE2, fill = "lightgrey", color = "lightgrey", size = 0.1) +
  geom_sf(data = globe_border_CAPTAIN2_EDGE2, fill = NA, color = "black", size = 0.5) +
  scale_fill_manual(values = c("Congruent high priority" = "red"),
                    name = "Priority areas") +
  labs(title = "Congruent high priority areas",
       subtitle = "Areas with priority > 0.9 in IUCN, FUSE and EDGE2 dimensions",
       x = NULL, y = NULL) +
  my_theme_CAPTAIN2_EDGE2

# Display the map
print(congruent_map)

ggsave(here::here("Outputs", "congruent_high_prioritiy_areas_09.png"), 
       congruent_map, 
       width = 10, height = 8, dpi = 300, bg = "white")

# Create detailed congruence categories
# Define logical conditions for each index
EDGE2_high <- combined_priorities$Priority_EDGE2 > 0.9
FUSE_high <- combined_priorities$Priority_FUSE > 0.9
IUCN_high <- combined_priorities$Priority_IUCN > 0.9

# Create detailed congruence categories
combined_priorities$congruence_category <- case_when(
  EDGE2_high & FUSE_high & IUCN_high ~ "All three indices",
  !EDGE2_high & FUSE_high & IUCN_high ~ "IUCN + FUSE",      # Changed from "FUSE + IUCN"
  EDGE2_high & !FUSE_high & IUCN_high ~ "IUCN + EDGE2",     # Changed from "EDGE2 + IUCN"
  EDGE2_high & FUSE_high & !IUCN_high ~ "FUSE + EDGE2",     # Changed from "EDGE2 + FUSE"
  !EDGE2_high & !FUSE_high & IUCN_high ~ "IUCN only",
  !EDGE2_high & FUSE_high & !IUCN_high ~ "FUSE only",
  EDGE2_high & !FUSE_high & !IUCN_high ~ "EDGE2 only",
  TRUE ~ "No high priority"
)

# Convert to factor with desired order
combined_priorities$congruence_category <- factor(
  combined_priorities$congruence_category,
  levels = c("All three indices",
             "IUCN + FUSE", "IUCN + EDGE2", "FUSE + EDGE2",  # Changed order
             "IUCN only", "FUSE only", "EDGE2 only")
)

# Filter data to only include high priority areas
high_priority_data <- combined_priorities[combined_priorities$congruence_category %in%
                                          c("All three indices", 
                                            "IUCN + FUSE", "IUCN + EDGE2", "FUSE + EDGE2",  # Changed names
                                            "IUCN only", "FUSE only", "EDGE2 only"), ]

# Debug: Check the data
cat("Data check:\n")

Data check:

cat("Total high priority areas:", nrow(high_priority_data), "\n")

Total high priority areas: 3139 

print(table(high_priority_data$congruence_category))

All three indices       IUCN + FUSE      IUCN + EDGE2      FUSE + EDGE2 
             1459                75              1130               275 
        IUCN only         FUSE only        EDGE2 only 
               29                17               154 

# Check if the factor levels are properly set
cat("\nFactor levels:\n")

Factor levels:

print(levels(high_priority_data$congruence_category))

[1] "All three indices" "IUCN + FUSE"       "IUCN + EDGE2"     
[4] "FUSE + EDGE2"      "IUCN only"         "FUSE only"        
[7] "EDGE2 only"       

# Check for any issues with the data
cat("\nData structure check:\n")

Data structure check:

print(str(high_priority_data$congruence_category))

 Factor w/ 7 levels "All three indices",..: 6 6 6 6 6 6 6 6 4 4 ...
NULL

# Test a simple plot first
cat("\nTesting simple plot...\n")

Testing simple plot...

test_plot <- ggplot() +
  geom_sf(data = high_priority_data, aes(fill = congruence_category)) +
  scale_fill_manual(values = c("All three indices" = "#8B0000",
                              "EDGE2 + FUSE" = "#FF4500",
                              "EDGE2 + IUCN" = "#FF6347",
                              "FUSE + IUCN" = "#FFA500",
                              "EDGE2 only" = "#4169E1",
                              "FUSE only" = "#32CD32",
                              "IUCN only" = "#9370DB")) +
  theme_void()

# Try to print the simple test plot
tryCatch({
  print(test_plot)
  cat("Simple plot works!\n")
}, error = function(e) {
  cat("Simple plot failed with error:", e$message, "\n")
})

Simple plot works!

# If simple plot fails, let's check the data more thoroughly
if (nrow(high_priority_data) == 0) {
  cat("No high priority data found! Checking original data...\n")
  cat("EDGE2 > 0.9:", sum(combined_priorities$Priority_EDGE2 > 0.9, na.rm = TRUE), "\n")
  cat("FUSE > 0.9:", sum(combined_priorities$Priority_FUSE > 0.9, na.rm = TRUE), "\n") 
  cat("IUCN > 0.9:", sum(combined_priorities$Priority_IUCN > 0.9, na.rm = TRUE), "\n")
}

# Extract coordinates for the congruence plot
if (nrow(high_priority_data) > 0) {
  coords <- st_coordinates(st_centroid(high_priority_data))
  plot_data <- data.frame(
    x = coords[,1],
    y = coords[,2], 
    category = high_priority_data$congruence_category
  )
  
  # Remove any rows with NA category
  plot_data <- plot_data[!is.na(plot_data$category), ]
  
  cat("Creating enhanced congruence plot...\n")
  cat("Plot data dimensions:", nrow(plot_data), "points\n")
  
  # Create the enhanced congruence plot
  congruence_map <- ggplot(plot_data, aes(x = x, y = y, color = category)) +
    geom_point(size = 1.2, alpha = 0.85, stroke = 0) +  # Larger points, no stroke to reduce overlap
    scale_color_manual(
      values = c("All three indices" = "#8B0000",      # Dark red
                "EDGE2 + FUSE" = "#FF8C00",            # Dark orange  
                "EDGE2 + IUCN" = "#DC143C",            # Crimson
                "FUSE + IUCN" = "#FFD700",             # Gold
                "EDGE2 only" = "#4169E1",              # Royal blue
                "FUSE only" = "#32CD32",               # Lime green
                "IUCN only" = "#9370DB"),              # Medium purple
      name = "High Priority\nCongruence",
      guide = guide_legend(
        override.aes = list(size = 4, alpha = 1),  # Larger, more opaque legend points
        title.position = "top",
        title.hjust = 0.5,
        ncol = 4  # Two columns for legend to save space
      )) +
    labs(
      title = "Global Conservation Priority Congruence Analysis",
      subtitle = "High priority areas (>0.9) showing agreement patterns across EDGE2, FUSE, and IUCN indices",
      x = "Longitude", 
      y = "Latitude",
      caption = paste("Total high priority areas:", nrow(plot_data))
    ) +
    theme_minimal() +
    theme(
      # Plot aesthetics
      plot.background = element_rect(fill = "white", color = NA),
      panel.background = element_rect(fill = "#f8f9fa", color = NA),
      panel.grid.major = element_line(color = "white", size = 0.5, linetype = "solid"),
      panel.grid.minor = element_line(color = "white", size = 0.25, linetype = "solid"),
      
      # Text styling
      plot.title = element_text(size = 16, hjust = 0.5, face = "bold", 
                               margin = margin(b = 10)),
      plot.subtitle = element_text(size = 12, hjust = 0.5, color = "gray30",
                                  margin = margin(b = 20)),
      plot.caption = element_text(size = 10, color = "gray50", hjust = 1),
      
      # Axis styling
      axis.title = element_text(size = 12, face = "bold"),
      axis.text = element_text(size = 10, color = "gray30"),
      axis.ticks = element_line(color = "gray50", size = 0.5),
      
      # Legend styling
      legend.position = "bottom",
      legend.background = element_rect(fill = "white", color = "gray80", size = 0.5),
      legend.margin = margin(t = 15, r = 10, b = 10, l = 10),
      legend.title = element_text(size = 12, face = "bold"),
      legend.text = element_text(size = 10),
      legend.key.size = unit(0.8, "cm"),
      
      # Panel border
      panel.border = element_rect(color = "gray70", fill = NA, size = 0.5)
    ) +
    # Add coordinate system and aspect ratio
    coord_fixed(ratio = 1.3) +  # Adjust ratio for better map appearance
    # Add subtle borders around plot area
    scale_x_continuous(expand = c(0.02, 0.02)) +
    scale_y_continuous(expand = c(0.02, 0.02))
  
  # Print the enhanced plot
  print(congruence_map)
  cat("Enhanced congruence plot created successfully!\n")
  
  # Print summary by category
  cat("\nBreakdown by congruence category:\n")
  category_counts <- table(plot_data$category)
  category_percentages <- round(prop.table(category_counts) * 100, 1)
  
  for(i in 1:length(category_counts)) {
    cat(sprintf("%-20s: %4d points (%4.2f%%)\n", 
                names(category_counts)[i], 
                category_counts[i], 
                category_percentages[i]))
  }
  
} else {
  cat("No valid high priority data to plot\n")
}

Creating enhanced congruence plot...
Plot data dimensions: 3139 points

Enhanced congruence plot created successfully!

Breakdown by congruence category:
All three indices   : 1459 points (46.50%)
IUCN + FUSE         :   75 points (2.40%)
IUCN + EDGE2        : 1130 points (36.00%)
FUSE + EDGE2        :  275 points (8.80%)
IUCN only           :   29 points (0.90%)
FUSE only           :   17 points (0.50%)
EDGE2 only          :  154 points (4.90%)

# Print summary of congruence patterns
cat("\nCongruence Pattern Summary:\n")

Congruence Pattern Summary:

congruence_summary <- table(combined_priorities$congruence_category, useNA = "ifany")
print(congruence_summary)

All three indices       IUCN + FUSE      IUCN + EDGE2      FUSE + EDGE2 
             1459                75              1130               275 
        IUCN only         FUSE only        EDGE2 only              <NA> 
               29                17               154                 6 

# Calculate percentages
total_high_priority <- sum(congruence_summary[names(congruence_summary) != "No high priority"], na.rm = TRUE)
congruence_percentages <- round(congruence_summary / total_high_priority * 100, 2)
cat("\nPercentages of high priority areas:\n")

Percentages of high priority areas:

print(congruence_percentages[names(congruence_percentages) != "No high priority"])

All three indices       IUCN + FUSE      IUCN + EDGE2      FUSE + EDGE2 
            46.48              2.39             36.00              8.76 
        IUCN only         FUSE only        EDGE2 only              <NA> 
             0.92              0.54              4.91                   

# ms style map -----
# Step 1: Check original data
print("=== STEP 1: CHECKING ORIGINAL DATA ===")

[1] "=== STEP 1: CHECKING ORIGINAL DATA ==="

print(paste("Number of rows in plot_data:", nrow(plot_data)))

[1] "Number of rows in plot_data: 3139"

print("First few rows of plot_data:")

[1] "First few rows of plot_data:"

print(head(plot_data))

          x       y  category
1 -767764.9 5762063 FUSE only
2 -727356.2 5762063 FUSE only
3 -772256.2 5713698 FUSE only
4 -731611.1 5713698 FUSE only
5 -776677.2 5665099 FUSE only
6 -735799.5 5665099 FUSE only

print("Summary of coordinates:")

[1] "Summary of coordinates:"

print(summary(plot_data[c("x", "y")]))

       x                   y           
 Min.   :-10914637   Min.   :-5217532  
 1st Qu.: -3784462   1st Qu.:-1511094  
 Median :  8930102   Median :  848691  
 Mean   :  4702046   Mean   :  599936  
 3rd Qu.: 12115021   3rd Qu.: 2810253  
 Max.   : 16383151   Max.   : 5762063  

print("Unique categories:")

[1] "Unique categories:"

print(unique(plot_data$category))

[1] FUSE only         FUSE + EDGE2      All three indices EDGE2 only       
[5] IUCN + EDGE2      IUCN + FUSE       IUCN only        
7 Levels: All three indices IUCN + FUSE IUCN + EDGE2 ... EDGE2 only

# Step 2: Transform to sf object
print("=== STEP 2: CREATING SF OBJECT ===")

[1] "=== STEP 2: CREATING SF OBJECT ==="

congruence_sf <- st_as_sf(
  plot_data, 
  coords = c("x", "y"), 
  crs = mcbryde_thomas_2  # Assuming WGS84
)
print(paste("Number of sf features created:", nrow(congruence_sf)))

[1] "Number of sf features created: 3139"

print("Original CRS:")

[1] "Original CRS:"

print(st_crs(congruence_sf))

Coordinate Reference System:
  User input: +proj=mbt_s 
  wkt:
PROJCRS["unknown",
    BASEGEOGCRS["unknown",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6326]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8901]]],
    CONVERSION["unknown",
        METHOD["PROJ mbt_s"]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]

print("Bounding box before projection:")

[1] "Bounding box before projection:"

print(st_bbox(congruence_sf))

     xmin      ymin      xmax      ymax 
-10914637  -5217532  16383151   5762063 

# Step 3: Project to McBryde-Thomas 2
print("=== STEP 3: PROJECTING TO MCBRYDE-THOMAS 2 ===")

[1] "=== STEP 3: PROJECTING TO MCBRYDE-THOMAS 2 ==="

congruence_sf <- st_transform(congruence_sf, crs = mcbryde_thomas_2)
print("CRS after projection:")

[1] "CRS after projection:"

print(st_crs(congruence_sf))

Coordinate Reference System:
  User input: +proj=mbt_s 
  wkt:
PROJCRS["unknown",
    BASEGEOGCRS["unknown",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6326]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8901]]],
    CONVERSION["unknown",
        METHOD["PROJ mbt_s"]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]

print("Bounding box after projection:")

[1] "Bounding box after projection:"

print(st_bbox(congruence_sf))

     xmin      ymin      xmax      ymax 
-10914637  -5217532  16383151   5762063 

# Step 4: Check world projection
print("=== STEP 4: CHECKING WORLD PROJECTION ===")

[1] "=== STEP 4: CHECKING WORLD PROJECTION ==="

world_projected_congruence <- st_transform(world, crs = mcbryde_thomas_2)
print("World bounding box:")

[1] "World bounding box:"

print(st_bbox(world_projected_congruence))

     xmin      ymin      xmax      ymax 
-18030969  -8669729  18199360   8324207 

# Step 5: Create globe border
print("=== STEP 5: CREATING GLOBE BORDER ===")

[1] "=== STEP 5: CREATING GLOBE BORDER ==="

globe_bbox <- rbind(c(-180, -90), c(-180, 90), 
                    c(180, 90), c(180, -90), c(-180, -90))
globe_border_congruence <- st_polygon(list(globe_bbox)) %>%
  st_sfc(crs = 4326) %>%
  st_sf(data.frame(rgn = 'globe', geom = .)) %>%
  smoothr::densify(max_distance = 0.5) %>%
  st_transform(crs = mcbryde_thomas_2)
print("Globe border bounding box:")

[1] "Globe border bounding box:"

print(st_bbox(globe_border_congruence))

     xmin      ymin      xmax      ymax 
-18373133  -8669783  18373133   8669783 

# Step 6: Check intersection
print("=== STEP 6: CHECKING INTERSECTION ===")

[1] "=== STEP 6: CHECKING INTERSECTION ==="

points_in_globe <- st_intersects(congruence_sf, globe_border_congruence, sparse = FALSE)
print(paste("Points within globe boundary:", sum(points_in_globe)))

[1] "Points within globe boundary: 3139"

print(paste("Percentage of points in globe:", round(sum(points_in_globe)/nrow(congruence_sf)*100, 2), "%"))

[1] "Percentage of points in globe: 100 %"

# Step 7: Extract coordinates for verification
print("=== STEP 7: EXTRACTING COORDINATES ===")

[1] "=== STEP 7: EXTRACTING COORDINATES ==="

coords <- st_coordinates(congruence_sf)
print("Projected coordinate summary:")

[1] "Projected coordinate summary:"

print(summary(coords))

       X                   Y           
 Min.   :-10914637   Min.   :-5217532  
 1st Qu.: -3784462   1st Qu.:-1511094  
 Median :  8930102   Median :  848691  
 Mean   :  4702046   Mean   :  599936  
 3rd Qu.: 12115021   3rd Qu.: 2810253  
 Max.   : 16383151   Max.   : 5762063  

print("First few projected coordinates:")

[1] "First few projected coordinates:"

print(head(coords))

             X       Y
[1,] -767764.9 5762063
[2,] -727356.2 5762063
[3,] -772256.2 5713698
[4,] -731611.1 5713698
[5,] -776677.2 5665099
[6,] -735799.5 5665099

# Step 8: Create theme
my_theme_congruence <- theme_minimal() +
  theme(
    legend.position = "bottom",
    legend.direction = "horizontal",
    legend.box = "vertical",
    legend.margin = margin(t = 20, r = 0, b = 0, l = 0),
    legend.title = element_text(margin = margin(b = 10)),
    panel.background = element_rect(fill = "white", color = NA),
    plot.background = element_rect(fill = "white", color = NA),
    panel.grid = element_blank(),
    axis.text = element_blank(),      # Add this line
    axis.ticks = element_blank()      # Add this line
  )

# Step 9: Create plot with debugging
print("=== STEP 9: CREATING PLOT ===")

[1] "=== STEP 9: CREATING PLOT ==="

congruence_map <- ggplot() +
  geom_sf(data = congruence_sf, aes(color = category), size = 1.2, alpha = 0.85) +
  geom_sf(data = world_projected_congruence, fill = "lightgrey", color = "lightgrey", size = 0.1) +
  geom_sf(data = globe_border_congruence, fill = NA, color = "lightgrey", size = 0.5) +
  scale_color_manual(
  values = c("All three indices" = "#8B0000",
            "IUCN + FUSE" = "#FFD700",           # Changed from "FUSE + IUCN"
            "IUCN + EDGE2" = "#DC143C",          # Changed from "EDGE2 + IUCN"
            "FUSE + EDGE2" = "#FF8C00",          # Changed from "EDGE2 + FUSE"
            "IUCN only" = "#9370DB",
            "FUSE only" = "#32CD32",
            "EDGE2 only" = "#4169E1"),
  breaks = c("All three indices",
             "IUCN + FUSE", 
             "IUCN + EDGE2", 
             "FUSE + EDGE2",
             "IUCN only", 
             "FUSE only", 
             "EDGE2 only"),              # Finally: FUSE + EDGE2
    name = "High priority congruence",
    guide = guide_legend(
      override.aes = list(size = 4, alpha = 1),
      title.position = "top", 
      title.hjust = 0.5,
      ncol = 4
    )) +
  labs(
    title = "Global conservation priority congruence analysis",
    subtitle = "High priority areas (>0.9) showing agreement patterns across IUCN, FUSE and EDGE2 prioritisations",
    x = NULL, 
    y = NULL
  ) +
  my_theme_congruence

ggsave(here::here("outputs", "congruent_high_prioritiy_areas_09_all_indices.png"), 
       congruence_map, 
       width = 10, height = 8, dpi = 300, bg = "white")

print("=== STEP 10: DISPLAYING PLOT ===")

[1] "=== STEP 10: DISPLAYING PLOT ==="

print("Plot created successfully")

[1] "Plot created successfully"

# Alternative test plot with just points to see if they're visible
print("=== CREATING TEST PLOT WITH JUST POINTS ===")

[1] "=== CREATING TEST PLOT WITH JUST POINTS ==="

test_plot <- ggplot() +
  geom_sf(data = congruence_sf, aes(color = category), size = 2, alpha = 1) +
  scale_color_manual(
    values = c("All three indices" = "#8B0000",      
              "EDGE2 + FUSE" = "#FF8C00",            
              "EDGE2 + IUCN" = "#DC143C",            
              "FUSE + IUCN" = "#FFD700",             
              "EDGE2 only" = "#4169E1",              
              "FUSE only" = "#32CD32",               
              "IUCN only" = "#9370DB")) +
  labs(title = "Test Plot - Points Only")

print("Test plot created")

[1] "Test plot created"

# Display both plots
print(congruence_map)

print(test_plot)

Congruence tests

# Statistical congruence analysis following the methods section
# First, extract priority values for the comparison
priority_threshold <- 0.9

# Create binary matrices for high-priority areas (>0.9)
fuse_high <- CAPTAIN2_FUSE_data_nonzero$Priority > priority_threshold
edge2_high <- CAPTAIN2_EDGE2_data_nonzero$Priority > priority_threshold
iucn_high <- CAPTAIN2_IUCN_data_nonzero$Priority > priority_threshold

# Calculate observed overlaps
observed_fuse_edge2 <- sum(fuse_high & edge2_high)
observed_fuse_iucn <- sum(fuse_high & iucn_high)
observed_edge2_iucn <- sum(edge2_high & iucn_high)
observed_all_three <- sum(fuse_high & edge2_high & iucn_high)

# Calculate minimum high-priority cells for each comparison
min_fuse_edge2 <- min(sum(fuse_high), sum(edge2_high))
min_fuse_iucn <- min(sum(fuse_high), sum(iucn_high))
min_edge2_iucn <- min(sum(edge2_high), sum(iucn_high))

# Calculate percentage overlaps (conservative estimate)
overlap_fuse_edge2 <- (observed_fuse_edge2 / min_fuse_edge2) * 100
overlap_fuse_iucn <- (observed_fuse_iucn / min_fuse_iucn) * 100
overlap_edge2_iucn <- (observed_edge2_iucn / min_edge2_iucn) * 100

# Randomization test function
randomization_test <- function(index1_high, index2_high, n_permutations = 999) {
  observed_overlap <- sum(index1_high & index2_high)
  min_cells <- min(sum(index1_high), sum(index2_high))
  observed_percentage <- (observed_overlap / min_cells) * 100
  
  # Generate random overlaps
  random_overlaps <- replicate(n_permutations, {
    # Randomly shuffle the spatial distribution of index1
    shuffled_index1 <- sample(index1_high)
    random_overlap <- sum(shuffled_index1 & index2_high)
    (random_overlap / min_cells) * 100
  })
  
  # Calculate p-value
  if (observed_percentage > mean(random_overlaps)) {
    p_value <- sum(random_overlaps >= observed_percentage) / n_permutations
  } else {
    p_value <- sum(random_overlaps <= observed_percentage) / n_permutations
  }
  
  return(list(
    observed_percentage = observed_percentage,
    expected_percentage = mean(random_overlaps),
    p_value = p_value,
    random_overlaps = random_overlaps
  ))
}

# Perform randomization tests
print("=== STATISTICAL CONGRUENCE ANALYSIS ===")

[1] "=== STATISTICAL CONGRUENCE ANALYSIS ==="

cat("Threshold for high-priority areas:", priority_threshold, "\n\n")

Threshold for high-priority areas: 0.9 

# FUSE vs EDGE2
test_fuse_edge2 <- randomization_test(fuse_high, edge2_high)

# FUSE vs IUCN
test_fuse_iucn <- randomization_test(fuse_high, iucn_high)

# EDGE2 vs IUCN
test_edge2_iucn <- randomization_test(edge2_high, iucn_high)

# Create a summary table
congruence_results <- data.frame(
  Comparison = c("FUSE vs EDGE2", "FUSE vs IUCN", "EDGE2 vs IUCN"),
  Observed_Overlap_Percent = c(test_fuse_edge2$observed_percentage,
                               test_fuse_iucn$observed_percentage,
                               test_edge2_iucn$observed_percentage),
  Expected_Overlap_Percent = c(test_fuse_edge2$expected_percentage,
                               test_fuse_iucn$expected_percentage,
                               test_edge2_iucn$expected_percentage),
  P_Value = c(test_fuse_edge2$p_value,
              test_fuse_iucn$p_value,
              test_edge2_iucn$p_value),
  Significance = c(
    ifelse(test_fuse_edge2$p_value < 0.001, "***",
           ifelse(test_fuse_edge2$p_value < 0.01, "**",
                  ifelse(test_fuse_edge2$p_value < 0.05, "*", "ns"))),
    ifelse(test_fuse_iucn$p_value < 0.001, "***",
           ifelse(test_fuse_iucn$p_value < 0.01, "**",
                  ifelse(test_fuse_iucn$p_value < 0.05, "*", "ns"))),
    ifelse(test_edge2_iucn$p_value < 0.001, "***",
           ifelse(test_edge2_iucn$p_value < 0.01, "**",
                  ifelse(test_edge2_iucn$p_value < 0.05, "*", "ns")))
  )
)

print(congruence_results)

     Comparison Observed_Overlap_Percent Expected_Overlap_Percent     P_Value
1 FUSE vs EDGE2                 94.56324                 94.90466 0.062062062
2  FUSE vs IUCN                 74.34031                 75.63198 0.002002002
3 EDGE2 vs IUCN                 95.18378                 94.83053 0.009009009
  Significance
1           ns
2           **
3           **

Relationship with fishing effort : bivariate maps

# Load your CAPTAIN data
IUCN_data <- CAPTAIN2_IUCN_data_nonzero  # Assuming this exists from your previous code
FUSE_data <- CAPTAIN2_FUSE_data_nonzero  # Assuming this exists from your previous code
EDGE2_data <- CAPTAIN2_EDGE2_data_nonzero # Assuming this exists from your previous code

# Process and prepare the data function
process_bivariate_data <- function(priority_data, fishing_data, min_priority = 0) {
  # Filter for cells with priority > min_priority
  priority_data_filtered <- priority_data %>%
    filter(Priority > min_priority)
  
  # Prepare fishing data
  fishing_data <- fishing_data %>%
    rename(Longitude = lon_05deg, Latitude = lat_05deg, FishingHours = predicted_fishing_hours)
  
  # Join datasets
  combined_data <- priority_data_filtered %>%
    left_join(fishing_data, by = c("Longitude", "Latitude"))
  
  # Handle NAs in fishing hours (replace with 0)
  combined_data$FishingHours[is.na(combined_data$FishingHours)] <- 0
  
  # Normalize priorities to 0-1 range and log transform fishing hours
  combined_data <- combined_data %>%
    mutate(
      Priority_Norm = Priority, # Assuming already in 0-1 range
      # Log transform fishing hours to better handle skewed distribution
      FishingHours_Log = log1p(FishingHours),
      # Normalize log fishing hours
      FishingHours_Norm = (FishingHours_Log - min(FishingHours_Log, na.rm = TRUE)) /
                         (max(FishingHours_Log, na.rm = TRUE) - min(FishingHours_Log, na.rm = TRUE))
    )
  
  # Set projection
  mcbryde_thomas_2 <- "+proj=mbt_s"
  
  # Transform to sf object
  data_sf <- st_as_sf(combined_data, coords = c("Longitude", "Latitude"), crs = 4326) %>%
    st_transform(crs = mcbryde_thomas_2)
  
  return(data_sf)
}

# Process the three datasets with fishing data (only cells with Priority > 0)
iucn_fishing_sf <- process_bivariate_data(IUCN_data, aggregated_data, min_priority = 0)
fuse_fishing_sf <- process_bivariate_data(FUSE_data, aggregated_data, min_priority = 0)
edge2_fishing_sf <- process_bivariate_data(EDGE2_data, aggregated_data, min_priority = 0)

# Create color palette - using the PurpleOr scheme
create_bivariate_palette <- function() {
  map_pal_raw <- bi_pal(pal = 'PurpleOr', dim = 4, preview = FALSE)
  map_pal_mtx <- matrix(map_pal_raw, nrow = 4, ncol = 4)
  map_pal_mtx[3, ] <- colorspace::lighten(map_pal_mtx[3, ], .1)
  map_pal_mtx[2, ] <- colorspace::lighten(map_pal_mtx[2, ], .2)
  map_pal_mtx[1, ] <- colorspace::lighten(map_pal_mtx[1, ], .3)
  map_pal_mtx[ , 3] <- colorspace::lighten(map_pal_mtx[ , 3], .1)
  map_pal_mtx[ , 2] <- colorspace::lighten(map_pal_mtx[ , 2], .2)
  map_pal_mtx[ , 1] <- colorspace::lighten(map_pal_mtx[ , 1], .3)
  map_pal_mtx[1, 1] <- '#ffffee'
  map_pal <- as.vector(map_pal_mtx) %>% setNames(names(map_pal_raw))
  return(map_pal)
}

map_pal <- create_bivariate_palette()

# Color mapping function
get_color <- function(priority, fishing) {
  priority_class <- cut(priority, breaks = c(-Inf, 0.25, 0.5, 0.75, Inf), labels = 1:4)
  fishing_class <- cut(fishing, breaks = c(-Inf, 0.25, 0.5, 0.75, Inf), labels = 1:4)
  return(map_pal[(as.numeric(fishing_class)-1)*4 + as.numeric(priority_class)])
}

# Apply colors to the datasets
iucn_fishing_sf$bivariate_color <- mapply(get_color,
                                          iucn_fishing_sf$Priority_Norm,
                                          iucn_fishing_sf$FishingHours_Norm)

fuse_fishing_sf$bivariate_color <- mapply(get_color,
                                          fuse_fishing_sf$Priority_Norm,
                                          fuse_fishing_sf$FishingHours_Norm)

edge2_fishing_sf$bivariate_color <- mapply(get_color,
                                          edge2_fishing_sf$Priority_Norm,
                                          edge2_fishing_sf$FishingHours_Norm)

# Get world map and project
world <- ne_countries(scale = "medium", returnclass = "sf")
mcbryde_thomas_2 <- "+proj=mbt_s"
world_projected <- st_transform(world, crs = mcbryde_thomas_2)

# Create globe border
lon_points <- rep(c(-180, 180), each = 100)
lat_points <- c(seq(-90, 90, length.out = 100), seq(90, -90, length.out = 100))
globe_outline <- data.frame(lon = lon_points, lat = lat_points) %>%
  st_as_sf(coords = c("lon", "lat"), crs = 4326) %>%
  st_combine() %>%
  st_cast("POLYGON") %>%
  st_transform(crs = mcbryde_thomas_2)

# Function to create individual bivariate maps with labels and titles
create_bivariate_map <- function(data_sf, label, title) {
  ggplot() +
    geom_sf(data = data_sf, aes(color = bivariate_color), size = 0.1, alpha = 1) +
    geom_sf(data = world_projected, fill = "lightgray", color = "lightgray") +
    geom_sf(data = globe_outline, fill = NA, color = "grey70", size = 0.5) +
    scale_color_identity() +
    coord_sf() +
    theme_minimal() +
    labs(x = NULL, y = NULL, title = title) +  # Added title
    theme(panel.grid = element_blank(),
          axis.text = element_blank(),
          axis.ticks = element_blank(),
          plot.title = element_text(size = 14, face = "bold", hjust = 0.5),  # Title styling
          plot.margin = margin(5, 5, 5, 5, "mm")) +
    # Add smaller label in top-left corner
    annotate("text", x = -Inf, y = Inf, label = label,
             hjust = -0.2, vjust = 1.2, size = 4, fontface = "bold")  # Reduced from 8 to 4
}

# Create the three maps with labels and titles
iucn_plot <- create_bivariate_map(iucn_fishing_sf, "(A)", "IUCN")
fuse_plot <- create_bivariate_map(fuse_fishing_sf, "(B)", "FUSE")  
edge2_plot <- create_bivariate_map(edge2_fishing_sf, "(C)", "EDGE2")

# Create a larger legend for better visibility
larger_legend <- bi_legend(pal = map_pal, dim = 4,
                          xlab = 'Conservation\npriority',
                          ylab = 'Fishing\neffort',
                          size = 4) +  # Increased from 3 to 4 for bigger legend
  theme(
    axis.title = element_text(size = 10, face = "bold"),  # Reduced from 12 to 10
    axis.text = element_blank(),
    legend.text = element_text(size = 8),  # Reduced from 10 to 8
    legend.title = element_text(size = 10, face = "bold"),  # Reduced from 12 to 10
    plot.margin = margin(10, 10, 10, 10, "mm")
  )

# Create layout matrix: 3 rows, 2 columns
# Column 1: maps (wider)
# Column 2: empty, legend (middle row), empty
layout_matrix <- rbind(
  c(1, NA),      # Row 1: Map A, empty
  c(2, 4),       # Row 2: Map B, legend
  c(3, NA)       # Row 3: Map C, empty
)

# Arrange all plots together
combined_plot <- grid.arrange(
  iucn_plot,
  fuse_plot,
  edge2_plot,
  larger_legend,
  layout_matrix = layout_matrix,
  heights = c(1, 1, 1),      # Equal height for all three map rows
  widths = c(3, 1)           # Maps take 3/4 width, legend takes 1/4
)

#To display the plot in the quarto file 
grid.arrange(
  iucn_plot,
  fuse_plot,
  edge2_plot,
  larger_legend,
  layout_matrix = layout_matrix,
  heights = c(1, 1, 1),      # Equal height for all three map rows
  widths = c(3, 1)           # Maps take 3/4 width, legend takes 1/4
)

# Save with dimensions better suited for A4 portrait
ggsave(here::here("Outputs", "Fig.4_All_indices_fishing_bivariate_maps_vertical.png"),
       combined_plot, width = 8.3, height = 11, dpi = 300, bg = "white")  # A4 portrait dimensions